#------------------------------------------------------------------------------
# Import libraries
#------------------------------------------------------------------------------

rm(list=ls())
library(LalRUtils)
libreq(tidyverse, data.table, zoo, tictoc, 
       fst, fixest, PanelMatch, patchwork,
       rio, magrittr, janitor, did, panelView, 
       ggiplot, tictoc, binsreg, interflex)
set.seed(42)
theme_set(lal_plot_theme())

#------------------------------------------------------------------------------





#------------------------------------------------------------------------------
# Define Paths
#------------------------------------------------------------------------------

# R studio
setwd( dirname(rstudioapi::getActiveDocumentContext()$path) )
# R default : unccoment if you use default R
# setwd(getSrcDirectory(function(){})[1])

#------------------------------------------------------------------------------



#------------------------------------------------------------------------------
# Load Data
#------------------------------------------------------------------------------

vcf <- fread("vcf_data_complete.csv", sep = ",")
setnames( vcf, "d", "D")
vcf_data <- copy( vcf )

#------------------------------------------------------------------------------




#------------------------------------------------------------------------------
# Managing Data
#------------------------------------------------------------------------------

vcf = vcf_data[year>= 1990]
vcf[, pref_bin := ntile(cover_1990, 10)]
vcf[, time_fra := year - 2008]
vcf[, D_fra := sch * (year >= 2008)]

vcf_data <- copy( vcf )
vcf_data[, `:=`(id_f = as.factor(cellid),
                sty_f = as.factor(styear),
                year_f = as.factor(year),
                D_f    = as.factor(D)
)]
vcf_data[, distance := min_dist_to_mine * 110]
intsamp = vcf_data[distance <= 100]

mining_all = interflex(Y = "forest_index", 
                       D = "D_f", 
                       X = "distance",
                       FE = c("id_f", "sty_f"), 
                       cl = "id_f",
                       data = intsamp, 
                       theme.bw = T,
                       Xlabel = 'Distance', 
                       Dlabel = "PESA",
                       Ylabel = 'Forest Index (vcf)',
                       cutoffs = seq(0, 100, 10),
                       estimator = 'binning', CI = TRUE )


# Get coefficients
est_table <- mining_all$est.bin[[1]]
ps_coef <- est_table[, 2]
ps_se <- est_table[, 3]
beta1 <- c()
j <- 1
for ( i in 1:10){
  beta1[ j ] <- ps_coef[ i ]
  j <- j + 1
  
  beta1[ j ] <- ps_se[ i ]
  j <- j + 1
}

m0_ps <- feols( forest_index ~ 1| id_f + sty_f , 
                cluster = "id_f",
                data = intsamp )

# Change name of coefficients
treatmap =c(
  # GFC
  "def_ha"       = "Annual Deforestation in Hectares",
  "village"      = "Village",
  "village[t]"   = "Village + Village TT",
  # VCF
  "forest_index" = "Forest cover index",
  "green_index"  = "Non-Forest Vegetation",
  "built_index"  = "Bare Ground Indices",
  "cellid"       = "Pixel",
  "id_f"         = "Pixel",
  "cellid[t]"    = "pixel + pixel TT",
  # both
  "yr"           = "Year",
  "year"         = "Year",
  "styear"       = "State $\\times$ Year",
  "sty_f"        = "State $\\times$ Year",
  "D"            = "PESA $\\times$ Scheduled",
  "block"        = "Block",
  "blk"          = "Block"
)

fitstat_register("n_new", function(x) summary( x )$nobs , 
                 "\\# Observations")


intsamp[ , never_treated := max(D) == 0, cellid ]
controls_pre1 = intsamp[never_treated == 1 & year < first_pesa_exposure, 
                        mean(forest_index)]
treat_pre1    = intsamp[never_treated == 0 & year < first_pesa_exposure, 
                        mean(forest_index)]

desc = "Treatment effects on Annual Deforestation by Distance to Mines"
fn = "figure8_regs"; lab = "tab:figure8_regs";

etable( m0_ps,
        style.tex = style.tex(main = "base", 
                              depvar.title = "", 
                              model.title = "", 
                              var.title = "\\midrule", 
                              slopes.title = "\\midrule \\emph{Time Trends}", 
                              yesNo = c("$\\checkmark$", "")),
        signif.code = NA,
        fixef_sizes = T, 
        fixef_sizes.simplify = F,
        fitstat = c( "n_new" ), 
        tex = TRUE,
        dict = treatmap,
        label = lab,
        title = glue::glue("{desc}"),
        extralines = list(
          "\\midrule Distance (0-10] $\\times$ PESA $\\times$ Scheduled" = c( beta1[ 1 ] ) ,
          " "                              = c( beta1[ 2 ] ) ,
          "Distance (10-20] $\\times$ PESA $\\times$ Scheduled" = c( beta1[ 3 ] ) ,
          " "                              = c( beta1[ 4 ] ) ,
          "Distance (20-30] $\\times$ PESA $\\times$ Scheduled" = c( beta1[ 5 ] ),
          " "                              = c( beta1[ 6 ] ) ,
          "Distance (30-40] $\\times$ PESA $\\times$ Scheduled" = c( beta1[ 7]),
          " "                              = c( beta1[ 8]),
          "Distance (40-50] $\\times$ PESA $\\times$ Scheduled" = c( beta1[ 9]),
          " "                              = c( beta1[10]),
          "Distance (50-60] $\\times$ PESA $\\times$ Scheduled" = c( beta1[ 11]),
          " "                              = c( beta1[12]),
          "Distance (60-70] $\\times$ PESA $\\times$ Scheduled" = c( beta1[ 13]),
          " "                              = c( beta1[14]),
          "Distance (70-80] $\\times$ PESA $\\times$ Scheduled" = c( beta1[ 15]),
          " "                              = c( beta1[16]),
          "Distance (80-90] $\\times$ PESA $\\times$ Scheduled" = c( beta1[ 17]),
          " "                              = c( beta1[18]),
          "Distance (90-100] $\\times$ PESA $\\times$ Scheduled" = c( beta1[ 19]),
          " "                              = c( beta1[20]),
          " " = c( " ", "\\midrule \\emph{Summary Statistics}", " " ),
          "Mean Pre-Y (Non-Sch)" = c( controls_pre1[ 1 ] ),
          "Mean Pre-Y (Sch )"    = c( treat_pre1[ 1 ] ),
          "Dataset"     = c( "VCF" ),
          "Timespan"    = c( "1995-2017" ) 
        ),
        file = "main_figure8_tablePanelB.tex", 
        replace = TRUE
)

#------------------------------------------------------------------------------